Bunching of fiuxons by the Cherenkov radiation in Josephson multilayers 
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A single magnetic fluxon moving at a high velocity in a Josephson multilayer {e.g., high-temperature 
superconductor such as BSCCO) can emit electromagnetic waves (Cherenkov radiation), which leads 
to formation of novel stable dynamic states consisting of several bunched fiuxons. We find such 
bunched states in numerical simulation in the simplest cases of two and three coupled junctions. At 
a given driving current, several different bunched states are stable and move at velocities that are 
higher than corresponding single-fluxon velocity. These and some of the more complex higher-order 
bunched states and transitions between them are investigated in detail. 

PACS: 74.50.+r, 74.80.Dm, 4f.60.Bq, 



I. INTRODUCTION 

In the recent years, a great deal of attention was at- 
tracted to different kinds of solid state multilayered sys- 
tems, e.g., artificial Josephson and magnetic multilay- 
ers, high-temperature superconductors (HTS) and per- 
ovskites, to name just a few. Multilayers are attractive 
because it is often possible to multiply a physical effect 
achieved in one layer by N (and sometimes by A 2 ), where 
N is the number of layers. This can be exploited for fab- 
rication of novel solid-state devices. In addition, multi- 
layered solid state systems show a variety of new physical 
phenomena which result from the interaction between in- 
dividual layers. 

In this article we focus on Josephson multilayer, the 
simplest example of which is a stack consisting of just 
two long Josephson junctions (LJJs). The results of our 
consideration can be applied to intrinsically layered HTS 
materialsGl, since the Josephson-stackjripdel has proved 
to be appropriate for n these structuresa □. 

In earlier papersETu it was shown that, in some cases, 
a fluxon (Josephson vortex) moving in one of the layers 
of the stack may emit electromagnetic (plasma) waves by 
means of the Cherenkov mechanism. The fluxon together 
with its Cherenkov radiation have a profile of a traveling 
wave, 4>{x — ut), having an oscillating gradually decaying 
tail. Such a wave profile generates an effective potential 
for another fluxon which can be added into the system. 
If the second fluxon is trapped in one of the minima of 
this traveling potential, we can get a bunched state of two 
fiuxons. In such a state, two fiuxons can stably move at 
a small constant distance one from another, which is not 
possible otherwise. Fiuxons of the same polarity usually 



repel each other, even being located in different layers. 

Similar bunched states were already found in a discrete 
Josephson transmission linesa, as well as in long Joseph- 
son junctions with the so-called due to the surface 
impedance of the superconductort3~ll3. The dynamics of 
conventional LJJ is described by the sine- Gordon equa- 
tion which does not allow the fluxon to move faster than 
the Swihart velocity and, therefore, the Cherenkov radi- 
ation never appears. In both cases mentioned above (the 
discrete system or the system with the /3-term) , the per- 
turbation of the sine-Gordon equation results in a mod- 
ified dispersion relation for Josephson plasma waves and 
appearance of an oscillating tail. This tail, in turn, re- 
sults in an attractive interaction between fiuxons, i.e., 
bunching. Nevertheless, the mere presence of an oscillat- 
ing tail is not a sufficient condition for bunching. 

In this paper, we investigate the problem of fluxon 
bunching in a system of two and three inductively cou- 
pled junctions with a primary state [1|0] (one fluxon in 
the top junction and no fluxon in the bottom one) or 
[0 1 1 1 0] (a fluxon only in the middle junction of a 3-fold 
stack). We show that bunching is possible for some fluxon 
configurations and specific range of parameters of the sys- 
tem. In addition, it is found that the bunched states radi- 
ate less than single-fluxon states, and therefore can move 
with a higher velocity. Section || presents the results of 
numerical simulations, in section |l| we discuss the ob- 
tained results and a feasibility of experimental observa- 
tion of bunched states. We also derive a simple analytical 
expression which show the possibility of the existence of 
bunched states. Section fv\ concludes the work. 
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II. NUMERICAL SIMULATIONS 

The system of equations which describes the dynamics 
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where S (— 1 < 5 < 0) is a dimensionless coupling con- 
stant, J = j A /j B is the ratio of the critical currents, 
while a and 7 = j/j A are the damping coefficient and 
normalized bias current, respectively, that are assumed 
to be the same in both LJJs. ft is also assumed that 
other parameters of the junctions, such as the effective 
magnetic thicknesses and, capacitances, are the same. As 
has been shown earliertffl, the Cherenkov radiation in a 
two-fold stack may take place only if the fluxon is mov- 
ing in the junction with smaller j c . We suppose in the 
following that the fluxon moves in LJJ" 4 , that implies 
J < 1. 

In the case N = 3, we impose the symmetry condition 
(f) A = <j) C , which is natural when the fluxon moves in the 
middle layer, and, thus, we can rewrite equations from 
Ref. [ [l3| in the form 
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Note the factor 2 in the last term on the l.h.s. of Eq. (|J). 
In the case of three coupled LJJs, we assume J = 1, 
since for more than two coupled junctions the Cherenkov 
radiation can be, obtained for a uniform stack with equal 
critical currentsD 



A. Numerical technique 

The numerical procedure works as follows. For a given 
set of the LJJs parameters, we compute the current- 
voltage characteristic (IVC) of the system, i.e., V A ' B {^). 
To calculate the voltages V A ' B for fixed values of 7, we 
simulate the dynamics of the phases jb A ' B (xA) by solving 
Eqs. (@) and (§) for TV = 2 or Eqs. (g) and (g) for N = 3, 
using the periodic boundary conditions: 



b A ' B {x 



L) 
L) 



b A > B (x 



Q)- 
0), 



2-kN 



A.B, 



(5) 
(6) 



where N A ' B is the number of fluxons trapped in LJJ" 4 ' 5 . 
In order to simulate a quasi-infinite system, we have cho- 
sen annular geometry with the length (circumference) of 
the junction L — 100. 



To solve the differential equations, we use an explicit 
method [expressing cf) A ' B (t + At) as a function of <p A ' B (t) 
and (f> A ' B (t — At)], treating <f) xx with a five-point, (j)u with 
a three-point, and 4>t with a two-point symmetric finite- 
difference scheme. The spatial and time steps used for 
the simulations were 8x = 0.025, St = 0.00625. After the 
simulation of the phase dynamics for T = 10 time units, 
we calculate the average dc voltages V A ' B for this time 
interval as 
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The dc voltage at point x can be defined as average num- 
ber of fluxons (the flux) passed through the junction at 
this point. Since the average fluxon density is not sin- 
gular in any point of the junction (otherwise the energy 
will grow infinitely) , we conclude that average dc voltage 
is the same for any point x. Therefore, for faster con- 
vergence of our averaging procedure, we can additionally 
average the phases (f> A ' B in (0) over the length of the 
stack. 

After the values of V A,B were found as per Eq. (|7]), 
the evolution of the phases (f> A ' B (x, t) is simulated further 
during LIT time units, the dc voltages V A ' B are calcu- 
lated for this new time interval and compared with the 
previously calculated values. We repeat such iterations 
further, increasing the time interval by a factor 1.1 until 
the difference in dc voltages \V(l.l n+1 T)-V(l.l n T)\ ob- 
tained in two subsequent iterations becomes less than an 
accuracy 5V = 10~ 4 . The particular factor 1.1 was found 
to be quite optimal and to provide for fast convergence, 
as well as more efficient averaging of low harmonics on 
each subsequent step. Very small value of this factor, 
e.g., 1.01 (recall that only the values greater than 1 have 
meaning), may result in a very slow convergence in the 
case when <j){t) contains harmonics with the period > T. 
Large values of the factor, e.g., > 2, would consume a lot 
of CPU time already during the second or third iteration 
and, hence, are not good for practical use. 

Once the voltage averaging for current 7 is complete, 
the current 7 is increased by a small amount $7 = 0.005 
to calculate the voltages at the next point of the IVC. We 
use a distribution of the phases (and their derivatives) 
achieved in the previous point of the IVC as the initial 
distribution for the following point. 

Further description of the software used for simulation 
can be found in Ref. 03. 



B. Two coupled junctions 

For simulation we chose the following parameters of 
the system: S — —0.5 to be close to the limit of intrin- 
sically layered HTS, J = 0.5 to let the fluxon accelerate 
above the c_ and develop Cherenkov radiation tail. The 
velocity c_ is the smallest of Swihart velocities of the 
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system. It characterizes the propagation of the out-of- 
phase mode of Josephson plasma waves. The value of 
a = 0.04 is chosen somewhat higher than, e.g., in (Nb- 
Al-AlO^Ar-Nb stacks. This choice is dictated by the need 
to keep the quasi-infinite approximation valid and satisfy 
the condition aL ^> 1. Smaller a requires very large L 
and, therefore, unaffordably long simulation times. So, 
we made a compromise and chose the above a value. 

First, we simulated the IVC u(j) in the [1 1 0] state, by 
sweeping 7 from up to 1 and making snapshots of phase 
gradients at every point of the IVC. This IVC is shown 
in Fig. ||(a), and the snapshot of the phase gradient at 
7 = 0.3 is presented in Fig. |l|(b). As one can see, the 
Cherenkov radiation tail, which is present for u > c_, 
has a sequence of minima where the second fluxon may 
be trapped. 

1. [1 + 1|0] state 

In order to create a two-fluxon bunched state and check 
its stability, we used the following "solution-engineering" 
procedure. By taking a snapshot of the phase profiles 
4>a,b{x) at the bias value 70 = 0.3, we constructed an 
ansatz for the bunched solution in the form 

4>%%(x) = 4> A , B (x) + 4>aA* + As), (8) 

where Ax is chosen so that the center of the trailing 
fluxon is placed at one of the minima of the Cherenkov 
tail. For example, to trap the trailing fluxon in the first, 
second and third well, we used Ax = 0.9, Ax = 2.4 
and Ax = 3.9, respectively. The phase distribution (and 
derivatives), constructed in this way, were used as the 
initial condition for solving Eqs. (|l|) and (Q) numerically. 
As the system relaxed to the desired state [1 + 1|0], we 
further traced ^(7) curve, varying 70 down to and up 
to 1. 

We accomplished this procedure for a set of Ace values, 
trying to trap the second fluxon in every well. Fig. |l](c) 
shows that a stable, tightly bunched state of two flux- 
ons is indeed possible. Actually, all the [1 + 1 1 0] states 
obtained this way have been found to be stable, and we 
were able to trace their IVCs up and down, starting from 
the initial value of the bias current 7 = 0.3. For the case 
when the trailing fluxon is trapped in the first, second 
and third minima, such IVCs are shown in Fig. 0. 

The most interesting feature of these curves is that 
they correspond to the velocity of the bunched state 
that is higher than that of the [1 1 0] state, at the same 
value of the bias current. Comparing solutions shown 
in Figs. [l](b) and |l|(c), we see that the amplitude of the 
trailing tail is smaller for the bunched state. This circum- 
stance suggests the following explanation to the fact that 
the observed velocity is higher in the state [1 + 1 1 0] than 
in the singlc-fluxon one. Because the driving forces act- 
ing on two fluxons in the bunched and unbunched states 



are the same, the difference in their velocities can be at- 
tributed only to the difference in the friction forces. The 
friction force acting on the fluxon in one junction is 

/+00 
4>x4>t dx , (9) 
-00 

and the same holds for the other junction. By just look- 
ing at Fig. |l|(b) and (c) it is rather difficult to tell in 
which case the friction force is larger, but accurate cal- 
culations using Eq. ([)]) and profiles from Fig. [l](b) and (c) 
show that the friction force acting on two fluxons with 
the tails shown in Fig. |l|(b) is somewhat higher than that 
for Fig. 0(c). This result is not surprising if one recalls 
that, to create the bunched state, we have shifted the 
[1|0] state by about half of the tail oscillation period rel- 
ative to the other single-fluxon state. Due to this, the 
tails of the two fluxons add up out of phase and partly 
cancel each other, making the tail's amplitude behind the 
fluxon in the bunched state lower than that in the [1 1 0] 
state. 

From Fig. |2| it is seen that every bunched state exists 
in a certain range of values of the bias current. If the 
current is decreased below some threshold value, fluxons 
dissociate and start moving apart, so that the interac- 
tion between them becomes exponentially small. When 
the trailing fluxon sits in a minimum of the Cherenkov 
tail sufficiently far from the leading fluxon, the IVC cor- 
responding to this bunched state is almost undistinguish- 
able from that of the [ 1 1 0] state, as the two fluxons ap- 
proach the limit when they do not interact. We have 
found that IVCs for M > 3, where M is the potential 
well's number, is indeed almost identical to that of the 
[l|0L.state. In contrast to bunching of fluxons in discrete 
LJjB, the transitions from one bunched state to another 
with different M do not take place in our system. Thus, 
we can say that the current range of a bunched state with 
smaller M "eclipses" the bunched states with larger M. 

The profiles of solutions found for various values of 
the bias current are shown in Fig. ||. We notice that 
at the bottom of the step corresponding to the bunched 
state the radiation tail is much weaker and fluxons are 
bunched tighter. This is a direct consequence of the fact 
that at lower velocities the radiation wavelength and the 
distance between minima becomes smaller, and so does 
the distance between the two fluxons. At a low bias cur- 
rent, the radiation wavelength and, hence, width of the 
potential wells become very small and incommensurable 
with the fluxon's width. Therefore, the fluxon does not 
fit into the well and the bunched states virtually disap- 
pear. 

2. [1 1 1] state 

The initial condition for this state was constructed in 
a similar fashion to the [1 + 1 1 0] one, but now using a 
cross-sum of the shifted and unshiftcd solutions: 
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$$(a:)=&Mi(aO + &M(a: + AaO. (10) 

If for the [1 + 1|0] state, Ax was « (A - |)M, M 
1.2.... then in the [1|1] state we have to take Ax ea 
AM. We can also take M = i.e., Ax = 0, which 
corresponds to the degenerate case of the in-phase [1|1] 
state. The stability of this state was investigated-ip detail 
analytically by Gr0nbech- Jensen and co-authorgl3, and is 
outside the scope of this paper. 

Our efforts to create a bound state [1 1 1] using the phase 
in the form ( |l(i| ) with M = 1,2... have not lead to any 
stable configuration of bunched fluxons with Ax ^= 0. 

3. Higher-order states 

Looking at the phase gradient profiles shown in Fig. ||, 
one notes that these profiles are qualitatively very similar 
to the original profile of the soliton with a radiation tail 
behind it [see Fig. |l|(b)], with the only difference that 
there are two bunched solitons with a tail. So, we can 
try to construct two pairs of bunched fluxons moving 
together, i.e., get a [2 + 2|0] bunched state. As before, 
the trapping of the trailing pair is possible in one of the 
minima of the tail generated by the leading pair. To 
construct such a double-bunched state we employ the 
initial conditions obtained using Eq. (||) at the bias point 
7o = 0.3, using the steady phase distribution obtained 
for the [2|0] state at 70 = 0.3. The shift Air was chosen 
in such a way that a pair of fluxons fits into one of the 
minima of the tail. We note that in this case we needed 
to vary Ax a little bit before we have achieved trapping 
of the trailing pair in a desired well. 

Simulations show that the obtained [2 + 2|0] states are 
stable and demonstrate an even higher velocity of the 
whole four-fluxon aggregate. The corresponding IVCs 
and profiles are shown in Fig. |Z](a) and (b) , respectively. 
Note that at 7 < 0.22 the bunched state [2 + 2|0] splits 
first into [1 + I2 + 13 + I3 10] state (the subscripts denote 
the well's number M, counting from the previous fluxon), 
and at still lower bias current, 7 < 0.2, they split into 
two independent [1 + I2 1 0] and [1 + I5 1 0] states. This 
two states move with slightly different velocities and can 
collide with each other due to the periodic nature of the 
system. As a result of collisions, these states ultimately 
undergo a transformation into two independent [1 + I5 10] 
states. As the bias decreases below w 0.1, the velocity u 
becomes smaller than c_ and the Cherenkov radiation 
tails disappear. At this point, each of the [1 + I5 1 0] 
states smoothly transforms into two independent [1 1 0] 
states. The interaction between these states is exponen- 
tially small, with a characteristic length ~ 1 (or, Aj in 
physical units). We note that the interaction between 
kinks in the region u > c_, where they have tails, also 
decreases exponentially, but with a larger characteristic 
length ~ a -1 . 

The procedure of constructing higher-order bunched 
states can be performed using different states as "build- 



ing blocks" . In particular, we also tried to form the 
[2 + 1 1 0] bunched state. Note that if two different states 
are taken as building blocks, we need to match their ve- 
locities, and, hence, the wave lengths of the tail. Thus, we 
have to combine two states at the same velocity, rather 
than at the same bias current. Since different states have 
their own velocity ranges, it is not always possible. As 
an example, we have constructed a [2 + 1 1 0] state out 
of [2|0] state at 7 = 0.15 and [1|0] state at 7 = 0.45 
using an ansatz similar to (^). These states have ap- 
proximately the same velocity u « 0.95 (see to Fig. |^). 
The constructed state was simulated, starting from the 
points 7 = 0.3 and 7 = 0.35, tracing IVC up and down as 
before. Depending on the bias current the system ends 
up in different states, namely in the state [1 + li + I2 10] 
for 70 = 0.3, or in the state [1 + U + li|0] = [3|0] for 
70 = 0.35. The IVCs of both states are shown in Fig. |J. 
The profiles of the phase gradients are shown in Fig. 0(c). 

Our attempts to construct the states with a higher 
number of bunched fluxons, e.g., [4+4|0], have failed since 
four fluxons do not fit into one well. We have concluded 
that such states immediately get converted into one of 
the lower-order states. 



C. Three coupled junctions 

We have performed numerical simulation of Eqs. (|^) 
and (§, using the same technique as described in the 
previous section. Our intention here is to study the 
3-junction case in which the fluxon is put in the mid- 
dle junction ( [0 1 1 1 0] state). All other parameters were 
the same as in the case of the two-junction system, ex- 
cept for the ratio of the critical currents J, which was 
taken equal to one. This simplest choice is made be- 
cause in a system of N > 2 coupled identical junctions 
the Cherenkov radiation appears in a [0| . . . |0|1|0| . . . |0] 
state for u > c_ sa 0.765 (this pertains to S = —0.5). 

Fig. H shows the IVCs of the original state [0|1|0], 
as well as IVCs of the bunched state [0 1 1 + 1 1 0] for 
M = 1, 2, 3. The profiles of the phase gradients at 
points A through D are shown in Fig. || Qualitatively, 
the bunching in the 3-fold system takes place in a sim- 
ilar fashion as that in the 2-fold system. Nevertheless, 
we did not succeed in creating a stable fluxon configura- 
tion with M = 3, although the stable states with other 
M were obtained. We would like to mention, that when 
the second fluxon was put in the second minimum of the 
potential to get the state with M = 2, the state with 
M = 1 has been finally established as a result of relax- 
ation. The same behavior was observed when we put the 
fluxon initially in the third minimum, the system ended 
up in the state [1 + 1 2 10] . For M > 4, the behavior was 
as usual. We tried to vary Ax smoothly, so that the cen- 
ter of the trailing fluxon would correspond to different 
positions between the second and fourth well, but in this 
case we did not succeed to get [1 + I2 1 0] state. 
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Following the same way as for two coupled junctions, 
we tried to construct [0+ 1|1|0 + 1] states. As in the case 
N = 2, these states were found unstable for any M > 0, 
e.g., they would split into [0|1 + 1 2 |0] and [1| - 1 1 1] . The 
state [0|2 + 2|0] was not stable either for M = 1, 2, 3 and 
the bias currents 70 = 0.20, 0.30, 0.35. 

The state [0|2 + 1 1 0] = [0|3|0], constructed by combin- 
ing the solutions for the [0 1 1 1 0] and [0|2|0] states moving 
with equal velocities was found to be stable when start- 
ing at 7 = 0.25 and sweeping bias current up and down. 
The dependence 14(7) is shown in Fig. || One may note, 
that for the states [0|2|0] and [0 1 3 1 0] the dependence is 
not smooth. Indeed, for these states the Cherenkov radia- 
tion tail is so long (~ L), that our annular system cannot 
simulate an infinitely long system, resulting in Cherenkov 
resonances whidaJnevitably appear in the system with a 
finite perimeter&Ll. 



1. Since non-oscillating tails result only in repulsion 
between fluxons, while the oscillating tail leads to 
mutual trapping, the condition 

Re( Pl )<\ Pr \ , (13) 

can be imposed to secure bunching. Here pi is the 
root of Eq. ([[2]) which describes the left (oscillat- 
ing) tail of the leading (right) fluxon, and p r is the 
root of Eq. (|lj) which describes the right (non- 
oscillating) tail of the trailing (left) fluxon. 

2. The relativistically contracted fluxon must fit into 
the minimum of the tail, i.e., 

7T / XL 2 

IMP) > \IW- 1 ' (14) 



III. ANALYSIS AND DISCUSSION 

Because of the non-linear nature of the bunching prob- 
lem, it is hardly tractable analytically. Therefore, we here 
present an approach in which we analyze the asymptotic 
behavior of the fluxon's front and trailing tails in the 
linear approximation. This technique is similar to that 
employed in Ref. |[ We assume that, at distances which 
are large enough in comparison with the fluxon's size, the 
fluxon's profile is exponentially decaying, 

<j)(x, t) oc exp[p(x — ut)] , (11) 

where p is a complex number which can be found by 
substituting this expression into Eqs. (0) and As a 
result we arrive at an equation 

Y~gz — p 2 u 2 — 1 — apu _ 
-T%2 T3s2 ~ V 2 v? ~ 7 ~ a P u 

(12) 

In general, this yields a 4-th order algebraic equation 
which always has 4 roots. If we want to describe a soliton 
moving from left to right with a radiation tail behind it, 
we have to find the values p among the four roots which 
adequately describe the front and rear parts of the soli- 
ton. Because the front (right) part of the soliton is not 
oscillating, it is described by Eq. ( pT| ) with real p < 0. 
The rear (left) part of the soliton is the oscillating tail, 
consequently it should be described by Eq. (^) with com- 
plex p having Re(p) > 0, the period of oscillations being 
determined by the imaginary part of p. Analyzing the 
4-th order equation, we conclude that the two necessary 
types of the roots coexist only for u > c_ , which is quite 
an obvious result. 

To analyze the possibility of bunched state formation, 
we consider two fluxons situated at some distance from 
each other. We propose the following two conditions for 
the two fluxons to form a bunched state: 



where 7r/Im(p) is half of the wavelength of the tail- 
forming radiation (the well's width), and the ex- 
pression on the r.h.s. of Eq. ( |l4"| ) approximately 
corresponds to the contraction of the fluxon at the 
trans-Swihart velocities. Although our system is 
not Lorentz invariant, numerical simulations show 
that the fluxon indeed shrinks (not up to zero) 
when approaching the Swihart velocity c_ from 
both sides. 

Following this approach, we have found that the second 
condition (Q) is always satisfied. The first condition ( |l3| ) 
gives the following result. Bunching is possible at u > 
Ub > c_. The value of ut, can be calculated numerically 
and for S = -0.5, J = 0.5, a = 0.04 it is u b = 0.837. 
Looking at Fig. ||, we see that this velocity corresponds to 
the bias point where the [1 + 1 a/ 1 0] states cease to exist. 
Thus, our crude approximation reasonably predicts the 
velocity range where the bunching is possible. 

IV. CONCLUSION 

In this work we have shown by means of numerical 
simulations that: 

• The emission of the Cherenkov plasma waves by a 
fluxon moving with high velocity creates an effec- 
tive potential with many wells, where other fluxons 
can be trapped. This mechanism leads to bunching 
between fluxons of the same polarity. 

• We have proved numerically that in the system of 
two and three coupled junctions the bunched states 
for the fluxons in the same junction such as [1 + 1|0], 
[1 + 2|0], [2 + 2 1 0] , [0|1 + 1|0] are stable. The states 
with fluxons in different junctions like [1|0+ 1] and 
[0+l|l|0+l] are numerically found unstable (except 
for the degenerated case M = 0, when [1|1] is a 
simple in- phase state). 



5 



• Bunched fluxons propagate at a substantially 
higher velocity than the corresponding free ones at 
the same bias current, because of lower losses per 
fluxon. 

• When decreasing the bias current, transitions be- 
tween the bunched states with different separations 
between fluxons were not found. This behavior dif- 
fers from what is known for the bunched states in a 
discrete systemEI. In addition, a splitting of multi- 
fluxon states into the states with smaller numbers 
of bunched fluxons is observed. 
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FIG. 1. (a) The current-velocity characteristic 11(7) for 
the fluxon moving in the [1|0] state (from left to right), (b) 
The profiles of the phase gradients <p^' B {x) in the state [1|0] at 
7 = 0.3, corresponding to the bias point A shown in Fig. (a). 
The Cherenkov tail, present at u > c_ ~ 0.817, has a set 
of minima where the second fluxon can be trapped, (c) The 
profiles of 4>^' B (x) in the state [2|0] at the same value 7 = 0.3 
as (b). Two fluxons shown in Fig. (c) are almost undistin- 
guishable. 

FIG. 2. Current-velocity characteristics of different 
bunched states [2 1 0] : the second fluxon is trapped in the first 
minimum of the tail (state [1 + 1 1 1 0] ) , the second minimum 
(state [1 + I2 1 0] ) , and the third minimum (state [1 + I3 1 0] ) . 
The 7(11) curve for the [1|0] state is shown for comparison. 
The phase-gradient profiles corresponding to the bias points 
A through D are shown in Fig. H. 
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FIG. 3. The profiles of the phase gradients (j>£' B (x) in the 
[2|0] states at the bias points A through D marked in Fig. |^. 



FIG. 4. (a) Current-velocity characteristics of the 
bunched states [4|0], [3)0], and [2 + 1 1 0] . Phase profiles of 
[4|0] state and [3|0] state at 7 = 0.3 are shown in (b) and (c), 
respectively. 



FIG. 5. Current-velocity characteristics of the state 
[0 1 1 1 0] , bunched state [0|1 + 1a/|0] for three different cases, 
M = 1, 2, 3, and the state [0 1 3 1 0] . The profiles of the Joseph- 
son p_hase gradients at the points A through D are shown in 
Fig. 



FIG. 6. The profiles of the Josephson phase gradients 
rf>x' B (x) in [0 1 1 + 1a/|0] states at the points A through D 
marked in Fig. |B|. 
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